This document integrates analysis of DRIFTS DPT files from OPUSprocessing.Rmd and Leila_code_modified.Rmd.
This code is modified from code by Stephany Soledad Chacon, “Box> Salk Institute Project> Salk Data> FTIR_Intact_decomposition_pots> FTIR_analysis.Rmd”. Using data in the DPT format, we can extrapolate sample attributes from the file names, then normalize and visualize spectra across replicates, crops, and timepoints. The width of the line bounding the spectra represents the range of the data.
# Check if folder exists and list files
if (!dir.exists(dpt_folder)) {
stop("Directory does not exist: ", dpt_folder)
}
# List and sort files
cat("Files found in directory:\n")## Files found in directory:
## [1] "DRIFTS_pot001_13Cwheat_wk0_root_250725.0.dpt"
## [2] "DRIFTS_pot012_13Csoy_wk10_root_250827.0.dpt"
## [3] "DRIFTS_pot029_13Cwheat_wk10_root_250827.0.dpt"
## [4] "DRIFTS_pot029_13Cwheat_wk40_root_250827.3.dpt"
## [5] "DRIFTS_pot030_13Csoy_wk10_root_250725.0.dpt"
## [6] "DRIFTS_pot030_13Csoy_wk40_root_250827.0.dpt"
## [7] "DRIFTS_pot032_13Cwheat_wk0_root_250725.0.dpt"
## [8] "DRIFTS_pot047_13Crice_wk0_root_250725.0.dpt"
## [9] "DRIFTS_pot050_13Csoy_wk0_root_250827.0.dpt"
## [10] "DRIFTS_pot052_13CnoPlant_wk30_root_250919.0.dpt"
## [11] "DRIFTS_pot052_13CnoPlant_wk40_root_250919.0.dpt"
## [12] "DRIFTS_pot053_13Cwheat_wk10_root_250725.0.dpt"
## [13] "DRIFTS_pot053_13Cwheat_wk40_root_250919.0.dpt"
## [14] "DRIFTS_pot055_13Crice_wk10_root_250725.0.dpt"
## [15] "DRIFTS_pot055_13Crice_wk10_root_qmark_250725.0.dpt"
## [16] "DRIFTS_pot055_13Crice_wk40_root_250827.0.dpt"
## [17] "DRIFTS_pot056_12Csoy_wk10_root_250919.0.dpt"
## [18] "DRIFTS_pot062_13Csoy_wk0_root_250725.0.dpt"
## [19] "DRIFTS_pot079_13Crice_wk0_root_250725.0.dpt"
## [20] "DRIFTS_pot080_13Csoy_wk10_root_250725.0.dpt"
## [21] "DRIFTS_pot080_13Csoy_wk40_root_250919.0.dpt"
## [22] "DRIFTS_pot086_12Csoy_wk0_root_250919.0.dpt"
## [23] "DRIFTS_pot087_13Crice_wk10_root_250827.0.dpt"
## [24] "DRIFTS_pot087_13Crice_wk40_root_250919.0.dpt"
## [25] "DRIFTS_pot094_13Csoy_wk10_root_250725.0.dpt"
## [26] "DRIFTS_pot094_13Csoy_wk10_soil_250630.0.dpt"
## [27] "DRIFTS_pot094_13Csoy_wk40_root_250919.0.dpt"
## [28] "DRIFTS_pot103_13Crice_wk0_root_recNMR_250725.0.dpt"
## [29] "DRIFTS_pot103_13Crice_wk0_root_vial1_250725.0.dpt"
## [30] "DRIFTS_pot103_13Crice_wk0_root_vial2_250725.2.dpt"
## [31] "DRIFTS_pot103_13Crice_wk0_root_vial3_250725.0.dpt"
## [32] "DRIFTS_pot107_12Crice_wk0_root_really13_250725.0.dpt"
## [33] "DRIFTS_pot109_13Cwheat_wk40_root_250919.0.dpt"
## [34] "DRIFTS_pot119_13Crice_wk10_root_250827.0.dpt"
## [35] "DRIFTS_pot119_13Crice_wk40_root_250919.0.dpt"
DPTfiles <- sort(list.files(dpt_folder, pattern = ".dpt"))
cat("Number of .dpt files:", length(DPTfiles), "\n")## Number of .dpt files: 35
# Read files individually and combine
dpt_data_list <- list()
for (i in seq_along(DPTfiles)) {
file_path <- file.path(dpt_folder, DPTfiles[i])
# Read as lines and parse manually (DPT files are space-delimited)
lines <- readLines(file_path, warn = FALSE)
lines <- lines[lines != ""] # Remove empty lines
# Split each line by whitespace
split_lines <- strsplit(lines, "\\s+")
# Convert to data frame
temp_data <- tryCatch({
data.frame(
wavenumber = as.numeric(sapply(split_lines, `[`, 1)),
absorbance = as.numeric(sapply(split_lines, `[`, 2)),
stringsAsFactors = FALSE
)
}, error = function(e) NULL)
if (!is.null(temp_data)) {
# Add filename column
temp_data$filename <- DPTfiles[i]
# Store in list
dpt_data_list[[i]] <- temp_data
}
}
# Remove NULL entries (failed reads)
dpt_data_list <- dpt_data_list[!sapply(dpt_data_list, is.null)]
# Combine all data
dpt_data <- bind_rows(dpt_data_list)
# Remove any rows with NA values
dpt_data <- dpt_data[complete.cases(dpt_data), ]
# Extract metadata from filenames
dpt_data <- dpt_data %>%
mutate(
crop = case_when(
str_detect(filename, "noPlant") ~ "noPlant",
str_detect(filename, "wheat") ~ "wheat",
str_detect(filename, "rice") ~ "rice",
str_detect(filename, "soy") ~ "soy",
TRUE ~ NA_character_
),
timepoint = str_extract(filename, "wk\\d+"),
sample_type = case_when(
str_detect(filename, "root") ~ "root",
str_detect(filename, "soil") ~ "soil",
TRUE ~ NA_character_
),
# Extract ID (pot number) from filename
ID = case_when(
str_detect(filename, "pot[0-9]+") ~ str_extract(filename, "pot[0-9]+") %>%
str_remove("pot") %>%
str_pad(width = 3, pad = "0"),
# Alternative pattern for files without 'pot' prefix
str_detect(filename, "DRIFTS_[0-9]+") ~ str_extract(filename, "DRIFTS_[0-9]+") %>%
str_remove("DRIFTS_") %>%
str_pad(width = 3, pad = "0"),
TRUE ~ NA_character_
)
) %>%
filter(!is.na(sample_type)) %>% # Remove files without clear sample type
filter(!str_detect(filename, "KBr")) %>% # Remove KBr files
# Filter out anomalous spectra
filter(!filename %in% c("DRIFTS_pot103_13Crice_wk0_root_recNMR_250725.0.dpt",
"DRIFTS_pot103_13Crice_wk0_root_vial1_250725.0.dpt",
"DRIFTS_pot079_13Crice_wk0_root_250725.0.dpt",
"DRIFTS_pot047_13Crice_wk0_root_250725.0.dpt",
"DRIFTS_pot086_12Csoy_wk0_root_250919.0.dpt",
"DRIFTS_pot056_12Csoy_wk10_root_250919.0.dpt",
"DRIFTS_pot030_13Csoy_wk10_root_250725.0.dpt"))
# save as csv for record-keeping
write.csv(dpt_data, file.path(outputs_folder, "dpt_data.csv"), row.names = FALSE)
# Filter for root samples only
dpt_data_roots <- dpt_data %>%
filter(sample_type == "root") %>%
filter(!is.na(crop))
cat("Root samples by crop:\n")## Root samples by crop:
##
## noPlant rice soy wheat
## 5598 27990 22392 19593
## Root samples by timepoint:
##
## wk0 wk10 wk30 wk40
## 19593 25191 2799 27990
This code is modified from code by Stephany Soledad Chacon, “Box> Salk Institute Project> Salk Data> FTIR_Intact_decomposition_pots> FTIR_analysis.Rmd”. It creates a new column with normalized absorbance values for each sample, using baseline correction (subtracting the minimum) followed by normalization to the maximum absorbance value.
# Normalize absorbance data by sample (same approach as OPUSprocessing.Rmd)
dpt_data_roots <- dpt_data_roots %>%
group_by(filename) %>%
mutate(
# Baseline correction - subtract minimum
baseline_corrected = absorbance - min(absorbance, na.rm = TRUE),
# Normalize to maximum
normalized_absorbance = baseline_corrected / max(baseline_corrected, na.rm = TRUE)
) %>%
ungroup()These plots are of normalized absorbance values. The thickness of the opaque line represents the range of the data. Annotations indicate integration regions used for calculations of simple plant, complex plant, and microbial organic matter proportions.
| Simple plant (aliphatic): | 2976 - 2898 cm-1 + 2870 - 2839 cm-1 |
| Microbial: | 1660 - 1580 cm-1 |
| Complex plant (aromatic): | 1550 - 1500 cm-1 |
# Generic function to create ridge plots for any timepoint
create_timepoint_ridge_plot <- function(data, timepoint_week, title = NULL) {
# Filter data for specific timepoint
filtered_data <- data %>%
filter(timepoint == timepoint_week) %>%
mutate(crop = factor(crop, levels = c("noPlant", "wheat", "rice", "soy")))
# Create title if not provided
if (is.null(title)) {
week_num <- str_extract(timepoint_week, "\\d+")
title <- paste("Roots at week", week_num)
}
# Determine y-position for annotations based on data
max_y <- max(filtered_data$normalized_absorbance)*(length(unique(filtered_data$crop))+3)
# Create the plot
ridge_plot <- filtered_data %>%
ggplot(mapping = aes(x = wavenumber,
y = crop,
height = normalized_absorbance,
group = crop,
fill = crop,
color = crop)) +
geom_density_ridges(stat = "identity",
scale = 3,
alpha = 0.25) +
geom_vline(xintercept = 1500, linetype = "dotted", linewidth = 0.3) +
geom_vline(xintercept = 1550, linetype = "dotted", linewidth = 0.3) +
geom_vline(xintercept = 1580, linetype = "dashed", linewidth = 0.3) +
geom_vline(xintercept = 1660, linetype = "dashed", linewidth = 0.3) +
geom_vline(xintercept = 2839, linetype = "longdash", linewidth = 0.3) +
geom_vline(xintercept = 2870, linetype = "longdash", linewidth = 0.3) +
geom_vline(xintercept = 2898, linetype = "longdash", linewidth = 0.3) +
geom_vline(xintercept = 2976, linetype = "longdash", linewidth = 0.3) +
xlim(3777, 422) +
scale_y_discrete(expand = expansion(mult = c(0.0, 0.56))) + # Less space below, more above
theme_classic() +
scale_fill_manual(values = crop_colors) +
scale_color_manual(values = crop_colors) +
ggtitle(title) +
# spacer
annotate("text", x = 1900, y = max_y + 0.2, label = " ",
size = 3, hjust = 0.5) +
# Add simple plant region annotation
annotate("rect", xmin = 2976, xmax = 2898, ymin = -Inf, ymax = Inf,
alpha = 0.2, fill = "#e3c8b1") +
annotate("rect", xmin = 2870, xmax = 2839, ymin = -Inf, ymax = Inf,
alpha = 0.2, fill = "#e3c8b1") +
annotate("text", x = 2560, y = max_y, label = "Simple Plant",
size = 3, hjust = 0.5) +
# Add microbial region annotation
annotate("rect", xmin = 1660, xmax = 1580, ymin = -Inf, ymax = Inf,
alpha = 0.2, fill = "#a8cbad") +
annotate("text", x = 1870, y = max_y + 0.11, label = "Microbial",
size = 3, hjust = 0.5) +
# Add complex plant region annotation
annotate("rect", xmin = 1550, xmax = 1500, ymin = -Inf, ymax = Inf,
alpha = 0.2, fill = "#beb5d5") +
annotate("text", x = 1190, y = max_y + 0.11, label = "Complex Plant",
size = 3, hjust = 0.5)
return(list(plot = ridge_plot, data = filtered_data))
}
# Generic function to create sample count tables
create_sample_count_table <- function(filtered_data, timepoint_week) {
week_num <- str_extract(timepoint_week, "\\d+")
col_name <- paste("Week", week_num, "samples")
counts <- filtered_data %>%
select(filename, crop) %>%
distinct() %>%
count(crop) %>%
arrange(match(crop, c("noPlant", "wheat", "rice", "soy"))) %>%
mutate(!!col_name := paste(n, " ", crop)) %>%
select(all_of(col_name))
return(counts)
}
# Create plots and tables for each timepoint
timepoints <- sort(unique(dpt_data_roots$timepoint))
for (tp in timepoints) {
# Create plot and get filtered data
result <- create_timepoint_ridge_plot(dpt_data_roots, tp)
# Print the plot
print(result$plot)
# Create and print sample count table with results='asis'
counts_table <- create_sample_count_table(result$data, tp)
cat(knitr::kable(counts_table, format = "html", escape = FALSE))
cat("<br><br>") # Add HTML line breaks for spacing
}| Week 0 samples |
|---|
| 2 wheat |
| 3 rice |
| 2 soy |
| Week 10 samples |
|---|
| 2 wheat |
| 4 rice |
| 3 soy |
| Week 30 samples |
|---|
| 1 noPlant |
| Week 40 samples |
|---|
| 1 noPlant |
| 3 wheat |
| 3 rice |
| 3 soy |
# Generic function to create line plots for any timepoint
create_timepoint_line_plot <- function(data, timepoint_week, title = NULL) {
# Filter data for specific timepoint
filtered_data <- data %>%
filter(timepoint == timepoint_week) %>%
mutate(crop = factor(crop, levels = c("noPlant", "wheat", "rice", "soy")))
# Create title if not provided
if (is.null(title)) {
week_num <- str_extract(timepoint_week, "\\d+")
title <- paste("Roots at week", week_num)
}
# Create the plot - offset lines vertically by crop like ridge plots
# Only include crops that are actually present in the filtered data
all_crop_levels <- c("noPlant", "wheat", "rice", "soy")
present_crops <- intersect(all_crop_levels, unique(filtered_data$crop))
len <- length(present_crops)
# Scale the data to match ridge plot proportions
max_abs <- max(filtered_data$normalized_absorbance, na.rm = TRUE)
# Determine y-position for annotations based on data
text_y <- len * 3 + 1.75
# Create offsets only for crops that are present
crop_offsets <- setNames(seq(0, (len-1) * 3, 3), present_crops)
# Add offset to data
filtered_data_offset <- filtered_data %>%
mutate(crop_offset = crop_offsets[as.character(crop)],
normalized_absorbance_offset = (normalized_absorbance*4.444 + crop_offset))
line_plot <- filtered_data_offset %>%
ggplot(mapping = aes(x = wavenumber,
y = normalized_absorbance_offset,
group = filename,
color = crop)) +
xlim(3777, 422) +
geom_line(alpha = 0.7) +
scale_y_continuous(
breaks = crop_offsets,
labels = names(crop_offsets) # Reduce top space from 0.56 to 0.15
) +
theme_classic() +
scale_fill_manual(values = crop_colors) +
scale_color_manual(values = crop_colors) +
ggtitle(title) +
labs(x = "Wavenumber (cm⁻¹)",
y = "",
color = "Crop") +
# Keep your existing vertical lines and annotations
geom_vline(xintercept = 1500, linetype = "dotted", linewidth = 0.3) +
geom_vline(xintercept = 1550, linetype = "dotted", linewidth = 0.3) +
geom_vline(xintercept = 1580, linetype = "dashed", linewidth = 0.3) +
geom_vline(xintercept = 1660, linetype = "dashed", linewidth = 0.3) +
geom_vline(xintercept = 2839, linetype = "longdash", linewidth = 0.3) +
geom_vline(xintercept = 2870, linetype = "longdash", linewidth = 0.3) +
geom_vline(xintercept = 2898, linetype = "longdash", linewidth = 0.3) +
geom_vline(xintercept = 2976, linetype = "longdash", linewidth = 0.3) +
# Add region annotations - positioned to match ridge plot style
# spacer
annotate("text", x = 1900, y = 10.201,
label = " ", size = 3, hjust = 0.5) +
# Add simple plant region annotation
annotate("rect", xmin = 2976, xmax = 2898, ymin = -Inf, ymax = Inf,
alpha = 0.2, fill = "#e3c8b1") +
annotate("rect", xmin = 2870, xmax = 2839, ymin = -Inf, ymax = Inf,
alpha = 0.2, fill = "#e3c8b1") +
annotate("text", x = 2560, y = text_y,
label = "Simple Plant", size = 3, hjust = 0.5) +
# Add microbial region annotation
annotate("rect", xmin = 1660, xmax = 1580, ymin = -Inf, ymax = Inf,
alpha = 0.2, fill = "#a8cbad") +
annotate("text", x = 1870, y = text_y,
label = "Microbial", size = 3, hjust = 0.5) +
# Add complex plant region annotation
annotate("rect", xmin = 1550, xmax = 1500, ymin = -Inf, ymax = Inf,
alpha = 0.2, fill = "#beb5d5") +
annotate("text", x = 1200, y = text_y,
label = "Complex Plant", size = 3, hjust = 0.5)
return(list(plot = line_plot, data = filtered_data))
}
# Generic function to create sample count tables
create_sample_count_table <- function(filtered_data, timepoint_week) {
week_num <- str_extract(timepoint_week, "\\d+")
col_name <- paste("Week", week_num, "samples")
counts <- filtered_data %>%
select(filename, crop) %>%
distinct() %>%
count(crop) %>%
arrange(match(crop, c("noPlant", "wheat", "rice", "soy"))) %>%
mutate(!!col_name := paste(n, " ", crop)) %>%
select(all_of(col_name))
return(counts)
}
# Create plots and tables for each timepoint
timepoints <- sort(unique(dpt_data_roots$timepoint))
for (tp in timepoints) {
# Create plot and get filtered data
result <- create_timepoint_line_plot(dpt_data_roots, tp)
# Print the plot
print(result$plot)
# Create and print sample count table with results='asis'
counts_table <- create_sample_count_table(result$data, tp)
cat(knitr::kable(counts_table, format = "html", escape = FALSE))
cat("<br><br>") # Add HTML line breaks for spacing
}| Week 0 samples |
|---|
| 2 wheat |
| 3 rice |
| 2 soy |
| Week 10 samples |
|---|
| 2 wheat |
| 4 rice |
| 3 soy |
| Week 30 samples |
|---|
| 1 noPlant |
| Week 40 samples |
|---|
| 1 noPlant |
| 3 wheat |
| 3 rice |
| 3 soy |
These plots are of normalized absorbance values. The thickness of the opaque line represents the range of the data. Annotations indicate integration regions used for calculations of simple plant, complex plant, and microbial organic matter proportions.
| Simple plant (aliphatic): | 2976 - 2898 cm-1 + 2870 - 2839 cm-1 |
| Microbial: | 1660 - 1580 cm-1 |
| Complex plant (aromatic): | 1550 - 1500 cm-1 |
# Define crop-specific color palettes
crop_timepoint_colors <- list(
wheat = c("wk0" = "#B4D586", "wk10" = "#98C559", "wk20" = "#79A63A",
"wk30" = "#58792A", "wk40" = "#374B1B"),
rice = c("wk0" = "#FDB35E", "wk10" = "#FC9D33", "wk29" = "#F18204",
"wk40" = "#B56203"),
soy = c("wk0" = "#FE71B1", "wk10" = "#FE318E", "wk20" = "#F4016E",
"wk30" = "#B70153", "wk40" = "#7A0137")
)
# Generic function to create crop-specific DRIFTS plots
create_crop_plot <- function(data, crop_name, color_palette) {
# Filter data for specific crop
crop_data <- data %>%
group_by(timepoint) %>%
filter(crop == crop_name)
# Ensure absorbance is numeric
crop_data$normalized_absorbance <- as.numeric(crop_data$normalized_absorbance)
# Create the plot
crop_plot <- ggplot(crop_data,
aes(x = wavenumber,
y = normalized_absorbance,
color = timepoint)) +
xlim(3777, 422) +
geom_line() +
theme_classic() +
scale_color_manual(values = color_palette) +
ggtitle(paste(str_to_title(crop_name), "Roots DRIFTS")) +
geom_vline(xintercept = 1500, linetype = "dotted", linewidth = 0.3) +
geom_vline(xintercept = 1550, linetype = "dotted", linewidth = 0.3) +
geom_vline(xintercept = 1580, linetype = "dashed", linewidth = 0.3) +
geom_vline(xintercept = 1660, linetype = "dashed", linewidth = 0.3) +
geom_vline(xintercept = 2839, linetype = "longdash", linewidth = 0.3) +
geom_vline(xintercept = 2870, linetype = "longdash", linewidth = 0.3) +
geom_vline(xintercept = 2898, linetype = "longdash", linewidth = 0.3) +
geom_vline(xintercept = 2976, linetype = "longdash", linewidth = 0.3) +
# Add simple plant region annotation
annotate("rect", xmin = 2976, xmax = 2898, ymin = -Inf, ymax = Inf,
alpha = 0.2, fill = "#e3c8b1") +
annotate("rect", xmin = 2870, xmax = 2839, ymin = -Inf, ymax = Inf,
alpha = 0.2, fill = "#e3c8b1") +
annotate("text", x = 2550, y = max(crop_data$normalized_absorbance), label = "Simple Plant",
size = 3, hjust = 0.5) +
# Add microbial region annotation
annotate("rect", xmin = 1660, xmax = 1580, ymin = -Inf, ymax = Inf,
alpha = 0.2, fill = "#a8cbad") +
annotate("text", x = 1850, y = max(crop_data$normalized_absorbance), label = "Microbial",
size = 3, hjust = 0.5) +
# Add complex plant region annotation
annotate("rect", xmin = 1550, xmax = 1500, ymin = -Inf, ymax = Inf,
alpha = 0.2, fill = "#beb5d5") +
annotate("text", x = 1200, y = max(crop_data$normalized_absorbance)+0.02, label = "Complex Plant",
size = 3, hjust = 0.5)
return(list(plot = crop_plot, data = crop_data))
}
# Generic function to create crop ID tables
create_crop_id_table <- function(crop_data, crop_name) {
combinations <- crop_data %>%
select(ID, timepoint) %>%
distinct() %>%
arrange(ID, timepoint) %>%
group_by(timepoint) %>%
summarise(pot_ids = {
ids <- paste(ID, sep = "")
# Group every 3 items and join with line breaks
grouped <- split(ids, ceiling(seq_along(ids)/5))
lines <- sapply(grouped, function(x) paste(x, collapse = ", "))
paste(lines, collapse = "<br>")
}, .groups = 'drop') %>%
pivot_wider(names_from = timepoint, values_from = pot_ids)
return(combinations)
}
# Generate plots and tables for each crop
crops <- c("wheat", "rice", "soy")
for (crop in crops) {
# Create plot and get data
result <- create_crop_plot(dpt_data_roots, crop, crop_timepoint_colors[[crop]])
# Print the plot
print(result$plot)
cat("<br>") # Add space between plot and table
# Create and print ID combinations table
id_table <- create_crop_id_table(result$data, crop)
cat(paste0("<h4>", str_to_title(crop), " pot root samples</h4>"))
cat(knitr::kable(id_table, format = "html", escape = FALSE, table.attr = 'style="border-collapse: separate; border-spacing: 35px 0;"'))
cat("<br><br>") # Add HTML line breaks for spacing
}| wk0 | wk10 | wk40 |
|---|---|---|
| 001, 032 | 029, 053 | 029, 053, 109 |
| wk0 | wk10 | wk40 |
|---|---|---|
| 103, 107 | 055, 087, 119 | 055, 087, 119 |
| wk0 | wk10 | wk40 |
|---|---|---|
| 050, 062 | 012, 080, 094 | 030, 080, 094 |
Modified from “Berhe_Ghezzehei_Lab/DRIFTS Data Analysis/Visualizing_FTIR_w_Elemental_Data.R” which starts by calculating aliphatic and aromatic indices, then deriving simple and complex plant and microbial proportions before visualizing everything. This procedure sums non-normalized absorbance values over a specified range, and divides by the total peak area in the ranges of interest to get the proportions.
# Function to perform spectral integration on DPT data (based on Leila_code_modified.Rmd)
perform_spectral_integration <- function(dpt_data) {
cat("Starting spectral integration processing...\n")
# Load required library for integration
if (!require(pracma, quietly = TRUE)) {
install.packages("pracma")
library(pracma)
}
# Prepare the data in the format needed for integration
# Rename columns to match expected format
comb <- dpt_data %>%
select(filename, wavenumber, absorbance) %>%
rename(sample = filename, wavelength = wavenumber, abs = absorbance)
# Define spectral windows of interest (same as Leila's code)
red1 <- comb[comb$wavelength > 2839 & comb$wavelength < 2870, ] # Aliphatic 1
red2 <- comb[comb$wavelength > 2898 & comb$wavelength < 2976, ] # Aliphatic 2
red3 <- comb[comb$wavelength > 1580 & comb$wavelength < 1660, ] # Microbial
red4 <- comb[comb$wavelength > 1500 & comb$wavelength < 1550, ] # Complex plant
# Get unique samples
sample <- unique(comb$sample)
# Create integration results data frame
ints <- data.frame(
sample = sample,
int_2839_2870 = numeric(length(sample)),
int_2898_2976 = numeric(length(sample)),
int_1580_1660 = numeric(length(sample)),
int_1500_1550 = numeric(length(sample)),
stringsAsFactors = FALSE
)
# Calculate integrated areas for each spectral window and sample
cat("Processing", length(sample), "samples for integration...\n")
for (i in seq_len(nrow(ints))) {
current_sample <- sample[i]
# Extract data for current sample from each spectral window
sample_red1 <- red1[red1$sample == current_sample, ]
sample_red2 <- red2[red2$sample == current_sample, ]
sample_red3 <- red3[red3$sample == current_sample, ]
sample_red4 <- red4[red4$sample == current_sample, ]
# Calculate area under curve (integration) using trapezoidal rule
if (nrow(sample_red1) > 1) {
ints$int_2839_2870[i] <- pracma::trapz(sample_red1$wavelength, sample_red1$abs)
}
if (nrow(sample_red2) > 1) {
ints$int_2898_2976[i] <- pracma::trapz(sample_red2$wavelength, sample_red2$abs)
}
if (nrow(sample_red3) > 1) {
ints$int_1580_1660[i] <- pracma::trapz(sample_red3$wavelength, sample_red3$abs)
}
if (nrow(sample_red4) > 1) {
ints$int_1500_1550[i] <- pracma::trapz(sample_red4$wavelength, sample_red4$abs)
}
}
# Use existing metadata from the input data (no need to re-extract)
cat("Using existing metadata from DPT data...\n")
basic_meta <- dpt_data %>%
select(filename, crop, timepoint, sample_type, ID) %>%
distinct() %>%
rename(sample = filename, type = sample_type)
# Combine metadata with integration results
ftir_data <- merge(basic_meta, ints, by = "sample", all = TRUE)
# Rename integration columns to match expected names
names(ftir_data)[names(ftir_data) == "int_1580_1660"] <- "int_microbial"
names(ftir_data)[names(ftir_data) == "int_1500_1550"] <- "int_complex_plant"
# Calculate derived metrics (same as original)
ftir_data$aliphatic <- ftir_data$int_2839_2870 + ftir_data$int_2898_2976
ftir_data$total_peak_area <- ftir_data$aliphatic + ftir_data$int_microbial + ftir_data$int_complex_plant
ftir_data$simple_plant_prop <- ftir_data$aliphatic / ftir_data$total_peak_area
ftir_data$complex_plant_prop <- ftir_data$int_complex_plant / ftir_data$total_peak_area
ftir_data$microbial_prop <- ftir_data$int_microbial / ftir_data$total_peak_area
ftir_data$simple_microbial_prop <- ftir_data$simple_plant_prop / ftir_data$microbial_prop
ftir_data$complex_microbial_prop <- ftir_data$complex_plant_prop / ftir_data$microbial_prop
# Remove individual aliphatic integration columns (keep only derived metrics)
ftir_data <- ftir_data %>%
select(-int_2839_2870, -int_2898_2976)
# Save the processed data to the output directory
today <- format(Sys.Date(), "%Y-%m-%d")
output_file <- file.path(outputs_folder, paste0("processed_int_data_", today, ".csv"))
write.csv(ftir_data, file = output_file, row.names = FALSE)
cat("Saved processed integration data to:", output_file, "\n")
cat("Integration processing completed successfully.\n")
return(ftir_data)
}# Try to find the most recent processed integration data file
ftir_file <- find_most_recent_processed_file(integrated_ftir_data)## Found 10 processed integration files
## Using most recent file: processed_int_data_2025-10-23.csv
if (!is.null(ftir_file) && file.exists(ftir_file)) {
ftir_data <- read_csv(ftir_file, show_col_types = FALSE)
cat("Successfully loaded FTIR integration data from:", basename(ftir_file), "\n")
cat("Dimensions:", dim(ftir_data), "\n")
cat("Columns:", names(ftir_data), "\n")
} else {
cat("No processed integration data found. Processing raw DPT data to generate integration results...\n")
# Check if we have the required dpt_data_roots from previous chunk
if (!exists("dpt_data_roots") || nrow(dpt_data_roots) == 0) {
stop("No DPT data available for integration. Please ensure the DPT data import was successful.")
}
# Process raw DPT data using the integrated function
ftir_data <- perform_spectral_integration(dpt_data_roots)
# Verify the processing was successful
if (is.null(ftir_data) || nrow(ftir_data) == 0) {
stop("Failed to process DPT data for integration.")
}
cat("Successfully processed", nrow(ftir_data), "samples from DPT data\n")
cat("Generated columns:", names(ftir_data), "\n")
}## Successfully loaded FTIR integration data from: processed_int_data_2025-10-23.csv
## Dimensions: 32 17
## Columns: source ID isotope crop timepoint type notes run_date int_microbial int_complex_plant aliphatic total_peak_area simple_plant_prop complex_plant_prop microbial_prop simple_microbial_prop complex_microbial_prop
# Calculate summary statistics by group
crop_stats <- ftir_data %>%
mutate(crop = factor(crop, levels = c("noPlant", "wheat", "rice", "soy"))) %>%
group_by(crop) %>%
summarise(
Mean_simple_plant_prop = mean(simple_plant_prop, na.rm = TRUE),
SE_simple_plant_prop = sd(simple_plant_prop, na.rm = TRUE) / sqrt(n()),
Mean_complex_plant_prop = mean(complex_plant_prop, na.rm = TRUE),
SE_complex_plant_prop = sd(complex_plant_prop, na.rm = TRUE) / sqrt(n()),
Mean_microbial_prop = mean(microbial_prop, na.rm = TRUE),
SE_microbial_prop = sd(microbial_prop, na.rm = TRUE) / sqrt(n()),
.groups = 'drop'
)
timepoint_stats <- ftir_data %>%
group_by(timepoint) %>%
summarise(
Mean_simple_plant_prop = mean(simple_plant_prop, na.rm = TRUE),
SE_simple_plant_prop = sd(simple_plant_prop, na.rm = TRUE) / sqrt(n()),
Mean_complex_plant_prop = mean(complex_plant_prop, na.rm = TRUE),
SE_complex_plant_prop = sd(complex_plant_prop, na.rm = TRUE) / sqrt(n()),
Mean_microbial_prop = mean(microbial_prop, na.rm = TRUE),
SE_microbial_prop = sd(microbial_prop, na.rm = TRUE) / sqrt(n()),
.groups = 'drop'
)
crop_timepoint_stats <- ftir_data %>%
mutate(crop = factor(crop, levels = c("noPlant", "wheat", "rice", "soy"))) %>%
group_by(crop, timepoint) %>%
summarise(
Mean_simple_plant_prop = mean(simple_plant_prop, na.rm = TRUE),
SE_simple_plant_prop = sd(simple_plant_prop, na.rm = TRUE) / sqrt(n()),
Mean_complex_plant_prop = mean(complex_plant_prop, na.rm = TRUE),
SE_complex_plant_prop = sd(complex_plant_prop, na.rm = TRUE) / sqrt(n()),
Mean_microbial_prop = mean(microbial_prop, na.rm = TRUE),
SE_microbial_prop = sd(microbial_prop, na.rm = TRUE) / sqrt(n()),
.groups = 'drop'
)This code is modified from “Visualizing_FTIR_w_Elemental_Data.R”
# Recreate the EXACT plotting function from Leila_code_modified.Rmd
create_proportion_plots <- function(stats_data, group_col, title_prefix) {
# Choose colors based on grouping variable
colors_to_use <- if(group_col == "crop") crop_colors else timepoint_colors
p1 <- ggplot(data = stats_data, aes_string(x = group_col, y = "Mean_simple_plant_prop", fill = group_col)) +
geom_col(alpha = 0.8) +
geom_errorbar(aes_string(ymin = "Mean_simple_plant_prop - SE_simple_plant_prop",
ymax = "Mean_simple_plant_prop + SE_simple_plant_prop"),
width = 0.2) +
theme_classic() +
scale_fill_manual(values = colors_to_use) +
labs(title = paste0(title_prefix, "Simple Plant OM"),
x = str_to_title(group_col),
y = "Simple Plant OM Proportion") +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1),
plot.margin = margin(20, 20, 20, 20),
axis.title = element_text(size = 12),
plot.title = element_text(size = 10))
p2 <- ggplot(data = stats_data, aes_string(x = group_col, y = "Mean_complex_plant_prop", fill = group_col)) +
geom_col(alpha = 0.8) +
geom_errorbar(aes_string(ymin = "Mean_complex_plant_prop - SE_complex_plant_prop",
ymax = "Mean_complex_plant_prop + SE_complex_plant_prop"),
width = 0.2) +
theme_classic() +
scale_fill_manual(values = colors_to_use) +
labs(title = paste0(title_prefix, "Complex Plant OM"),
x = str_to_title(group_col),
y = "Complex Plant OM Proportion") +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1),
plot.margin = margin(20, 20, 20, 20),
axis.title = element_text(size = 12),
plot.title = element_text(size = 10))
p3 <- ggplot(data = stats_data, aes_string(x = group_col, y = "Mean_microbial_prop", fill = group_col)) +
geom_col(alpha = 0.8) +
geom_errorbar(aes_string(ymin = "Mean_microbial_prop - SE_microbial_prop",
ymax = "Mean_microbial_prop + SE_microbial_prop"),
width = 0.2) +
theme_classic() +
scale_fill_manual(values = colors_to_use) +
labs(title = paste0(title_prefix, "Microbial OM"),
x = str_to_title(group_col),
y = "Microbial OM Proportion") +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1),
plot.margin = margin(20, 20, 20, 20),
axis.title = element_text(size = 12),
plot.title = element_text(size = 10))
list(p1 = p1, p2 = p2, p3 = p3)
}
# Create EXACT plots by crop (from Leila_code_modified.Rmd)
crop_plots <- create_proportion_plots(crop_stats, "crop", "")
crop_combined <- plot_grid(crop_plots$p1, crop_plots$p2, crop_plots$p3, nrow = 1, labels = "auto")
print(crop_combined)# Create EXACT plots by timepoint (from Leila_code_modified.Rmd)
timepoint_plots <- create_proportion_plots(timepoint_stats, "timepoint", "")
timepoint_combined <- plot_grid(timepoint_plots$p1, timepoint_plots$p2, timepoint_plots$p3, nrow = 1, labels = "auto")
print(timepoint_combined)# Recreate the EXACT interaction plot (crop x timepoint) from Leila_code_modified.Rmd
interaction_plot <- ggplot(crop_timepoint_stats) +
geom_col(aes(x = crop, y = Mean_simple_plant_prop, fill = timepoint),
position = "dodge", alpha = 0.8) +
geom_errorbar(aes(x = crop,
ymin = Mean_simple_plant_prop - SE_simple_plant_prop,
ymax = Mean_simple_plant_prop + SE_simple_plant_prop,
group = timepoint),
position = position_dodge(0.9), width = 0.2) +
theme_classic() +
scale_fill_manual("Timepoint", values = timepoint_colors) +
labs(title = "Simple Plant OM by Crop and Timepoint",
x = "Crop Type",
y = "Simple Plant OM Proportion") +
theme(legend.position = "top",
plot.margin = margin(20, 20, 20, 20),
axis.text.x = element_text(angle = 45, hjust = 1))
print(interaction_plot)# Recreate EXACT stacked bar plots from Leila_code_modified.Rmd
create_stacked_data <- function(stats_data, group_col) {
# Reshape data for stacking (exact reproduction)
long_data <- stats_data %>%
select(all_of(group_col), Mean_simple_plant_prop, Mean_complex_plant_prop, Mean_microbial_prop) %>%
pivot_longer(cols = starts_with("Mean_"),
names_to = "func_type",
values_to = "proportion") %>%
mutate(func_type = case_when(
func_type == "Mean_simple_plant_prop" ~ "Simple Plant",
func_type == "Mean_complex_plant_prop" ~ "Complex Plant",
func_type == "Mean_microbial_prop" ~ "Microbial"
)) %>%
mutate(func_type = factor(func_type, levels = c("Complex Plant", "Simple Plant", "Microbial")))
long_data
}
# Create EXACT stacked plots
crop_stacked_data <- create_stacked_data(crop_stats, "crop")
crop_stacked_data$crop <- factor(crop_stacked_data$crop, levels = c("noPlant", "wheat", "rice", "soy"))
timepoint_stacked_data <- create_stacked_data(timepoint_stats, "timepoint")
p_crop_stacked <- ggplot(crop_stacked_data, aes(x = crop, y = proportion * 100, fill = func_type)) +
geom_col(position = "stack") +
scale_fill_manual(values = c("#93deed", "#d2f0f4", "#1ca5cf"),
breaks = c("Microbial", "Simple Plant", "Complex Plant")) +
theme_classic() +
labs(x = "Crop Type",
y = "Functional Group Amount (%)",
fill = "") +
theme(legend.position = "top",
plot.margin = margin(20, 27, 20, 13),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.margin = margin(10, 16, 10, 8))
p_timepoint_stacked <- ggplot(timepoint_stacked_data, aes(x = timepoint, y = proportion * 100, fill = func_type)) +
geom_col(position = "stack") +
scale_fill_manual(values = c("#e8aa9c", "#f5d8d1", "#c54e30"),
breaks = c("Microbial", "Simple Plant", "Complex Plant")) +
theme_classic() +
labs(x = "Timepoint",
y = "Functional Group Amount (%)",
fill = "") +
theme(legend.position = "top",
plot.margin = margin(20, 27, 20, 13),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.margin = margin(10, 16, 10, 8))
combined_stacked <- plot_grid(
p_crop_stacked,
p_timepoint_stacked,
nrow = 1,
labels = "auto"
)
print(combined_stacked)These ridge plots show the same spectral data as the original Ridge Plots section, but with annotations highlighting suberin-related peaks identified in the literature. The annotations indicate peaks associated with suberin content in plant roots, based on White et al. (2011, 2016).
| Aliphatic suberin moiety: | 2976, 2959, 2930, 2856 cm-1 (within 3000-2800 cm-1 range) |
| Suberin esters: | 1737-1740 cm-1 |
| Aromatic suberin moiety: | 1506 cm-1 (also lignin) |
White, K. E., Reeves, J. B., & Coale, F. J. (2011). Mid-infrared
diffuse reflectance spectroscopy for the rapid analysis of plant root
composition. Geoderma, 167–168, 197–203. https://doi.org/10.1016/j.geoderma.2011.08.009
or
White, K. E., Reeves, J. B., & Coale, F. J. (2016). Cell wall
compositional changes during incubation of plant roots measured by
mid-infrared diffuse reflectance spectroscopy and fiber analysis.
Geoderma, 264, 205–213. https://doi.org/10.1016/j.geoderma.2015.10.018
# Generic function to create ridge plots with suberin annotations
create_suberin_ridge_plot <- function(data, timepoint_week, title = NULL) {
# Filter data for specific timepoint
filtered_data <- data %>%
filter(timepoint == timepoint_week) %>%
mutate(crop = factor(crop, levels = c("noPlant", "wheat", "rice", "soy")))
# Create title if not provided
if (is.null(title)) {
week_num <- str_extract(timepoint_week, "\\d+")
title <- paste("Root suberin at week", week_num)
}
# Determine y-position for annotations based on data
max_y <- max(filtered_data$normalized_absorbance)*(length(unique(filtered_data$crop))+3)
# Create the plot
ridge_plot <- filtered_data %>%
ggplot(mapping = aes(x = wavenumber,
y = crop,
height = normalized_absorbance,
group = crop,
fill = crop,
color = crop)) +
geom_density_ridges(stat = "identity",
scale = 3,
alpha = 0.25) +
# Suberin peak annotations
# Aliphatic suberin region (3000-2800 cm-1)
annotate("rect", xmin = 2800, xmax = 3000, ymin = -Inf, ymax = Inf,
alpha = 0.1, fill = "#5076ab") +
# Specific aliphatic peaks
geom_vline(xintercept = 2976, linetype = "longdash", linewidth = 0.3, color = "#3D6AA9") +
geom_vline(xintercept = 2959, linetype = "longdash", linewidth = 0.3, color = "#3D6AA9") +
geom_vline(xintercept = 2930, linetype = "longdash", linewidth = 0.3, color = "#3D6AA9") +
geom_vline(xintercept = 2856, linetype = "longdash", linewidth = 0.3, color = "#3D6AA9") +
# Suberin esters (1737-1740 cm-1)
annotate("rect", xmin = 1733, xmax = 1741, ymin = -Inf, ymax = Inf,
alpha = 0.15, fill = "#4d653a") +
geom_vline(xintercept = 1737, linetype = "dashed", linewidth = 0.3, color = "#4d653a") +
# Aromatic suberin (1506 cm-1)
geom_vline(xintercept = 1506, linetype = "dotted", linewidth = 0.3, color = "#5f3e38") +
xlim(3200, 1200) +
scale_y_discrete(expand = expansion(mult = c(0.01, 0.56))) + # Less space below, more above
theme_classic() +
scale_fill_manual(values = crop_colors) +
scale_color_manual(values = crop_colors) +
ggtitle(title) +
# spacer
annotate("text", x = 1900, y = max_y + 0.2, label = " ",
size = 3, hjust = 0.5) +
# Add aliphatic annotation
annotate("text", x = 2600, y = max_y-0.2, label = "aliphatic suberin",
size = 3, hjust = 0.5) +
# Add esters annotation
annotate("text", x = 1910, y = max_y, label = "suberin esters",
size = 3, hjust = 0.5) +
# Add aromatic annotation
annotate("text", x = 1300, y = max_y + 0.1, label = "aromatic suberin",
size = 3, hjust = 0.5)
return(list(plot = ridge_plot, data = filtered_data))
}
# Get unique timepoints and sort them (use dpt_data_roots, not ftir_data)
timepoints <- sort(unique(dpt_data_roots$timepoint))
# Create plots and tables for each timepoint
timepoints <- sort(unique(dpt_data_roots$timepoint))
for (tp in timepoints) {
# Create plot and get filtered data
result <- create_suberin_ridge_plot(dpt_data_roots, tp)
# Print the plot
print(result$plot)
# Create and print sample count table with results='asis'
counts_table <- create_sample_count_table(result$data, tp)
cat(knitr::kable(counts_table, format = "html", escape = FALSE))
cat("<br><br>") # Add HTML line breaks for spacing
}| Week 0 samples |
|---|
| 2 wheat |
| 3 rice |
| 2 soy |
| Week 10 samples |
|---|
| 2 wheat |
| 4 rice |
| 3 soy |
| Week 30 samples |
|---|
| 1 noPlant |
| Week 40 samples |
|---|
| 1 noPlant |
| 3 wheat |
| 3 rice |
| 3 soy |
# Define crop-specific color palettes
crop_timepoint_colors <- list(
noPlant = c("wk0" = "#B68ADE", "wk10" = "#9C6BC7", "wk20" = "#7B41A9",
"wk30" = "#6f3e94", "wk40" = "#3F1C59"),
wheat = c("wk0" = "#B4D586", "wk10" = "#98C559", "wk20" = "#79A63A",
"wk30" = "#58792A", "wk40" = "#374B1B"),
rice = c("wk0" = "#FDB35E", "wk10" = "#FC9D33", "wk29" = "#F18204",
"wk40" = "#B56203"),
soy = c("wk0" = "#FE71B1", "wk10" = "#FE318E", "wk20" = "#F4016E",
"wk30" = "#B70153", "wk40" = "#7A0137")
)
# Generic function to create crop-specific DRIFTS plots with suberin annotations
create_crop_plot <- function(data, crop_name, color_palette) {
# Filter data for specific crop
crop_data <- data %>%
group_by(timepoint) %>%
filter(crop == crop_name)
# Ensure absorbance is numeric
crop_data$normalized_absorbance <- as.numeric(crop_data$normalized_absorbance)
# Create the plot
crop_plot <- ggplot(crop_data,
aes(x = wavenumber,
y = normalized_absorbance,
color = timepoint)) +
xlim(3200, 1200) +
geom_line() +
theme_classic() +
scale_color_manual(values = color_palette) +
ggtitle(paste(str_to_title(crop_name), "root suberin")) +
# Suberin peak annotations
# Aliphatic suberin region (3000-2800 cm-1)
annotate("rect", xmin = 2800, xmax = 3000, ymin = -Inf, ymax = Inf,
alpha = 0.1, fill = "#5076ab") +
# Specific aliphatic peaks
geom_vline(xintercept = 2976, linetype = "longdash", linewidth = 0.3, color = "#3D6AA9") +
geom_vline(xintercept = 2959, linetype = "longdash", linewidth = 0.3, color = "#3D6AA9") +
geom_vline(xintercept = 2930, linetype = "longdash", linewidth = 0.3, color = "#3D6AA9") +
geom_vline(xintercept = 2856, linetype = "longdash", linewidth = 0.3, color = "#3D6AA9") +
# Suberin esters (1737-1740 cm-1)
annotate("rect", xmin = 1733, xmax = 1741, ymin = -Inf, ymax = Inf,
alpha = 0.15, fill = "#4d653a") +
geom_vline(xintercept = 1737, linetype = "dashed", linewidth = 0.3, color = "#4d653a") +
# Aromatic suberin (1506 cm-1)
geom_vline(xintercept = 1506, linetype = "dotted", linewidth = 0.3, color = "#5f3e38") +
# Add aliphatic annotation
annotate("text", x = 2600, y = max(crop_data$normalized_absorbance)-0.05, label = "aliphatic suberin",
size = 3, hjust = 0.5) +
# Add esters annotation
annotate("text", x = 1910, y = max(crop_data$normalized_absorbance), label = "suberin esters",
size = 3, hjust = 0.5) +
# Add aromatic annotation
annotate("text", x = 1300, y = max(crop_data$normalized_absorbance)+0.02, label = "aromatic suberin",
size = 3, hjust = 0.5)
return(list(plot = crop_plot, data = crop_data))
}
# Generate plots and tables for each crop
crops <- c("noPlant", "wheat", "rice", "soy")
for (crop in crops) {
# Create plot and get data
result <- create_crop_plot(dpt_data_roots, crop, crop_timepoint_colors[[crop]])
# Print the plot
print(result$plot)
cat("<br>") # Add space between plot and table
# Create and print ID combinations table
id_table <- create_crop_id_table(result$data, crop)
cat(paste0("<h4>", str_to_title(crop), " pot root samples</h4>"))
cat(knitr::kable(id_table, format = "html", escape = FALSE, table.attr = 'style="border-collapse: separate; border-spacing: 35px 0;"'))
cat("<br><br>") # Add HTML line breaks for spacing
}| wk30 | wk40 |
|---|---|
| 052 | 052 |
| wk0 | wk10 | wk40 |
|---|---|---|
| 001, 032 | 029, 053 | 029, 053, 109 |
| wk0 | wk10 | wk40 |
|---|---|---|
| 103, 107 | 055, 087, 119 | 055, 087, 119 |
| wk0 | wk10 | wk40 |
|---|---|---|
| 050, 062 | 012, 080, 094 | 030, 080, 094 |
# Generic function to create line plots for any crop with suberin annotations
# Define crop-specific color palettes
crop_timepoint_colors <- list(
noPlant = c("wk0" = "#B68ADE", "wk10" = "#9C6BC7", "wk20" = "#7B41A9",
"wk30" = "#6f3e94", "wk40" = "#3F1C59"),
wheat = c("wk0" = "#B4D586", "wk10" = "#98C559", "wk20" = "#79A63A",
"wk30" = "#58792A", "wk40" = "#374B1B"),
rice = c("wk0" = "#FDB35E", "wk10" = "#FC9D33", "wk29" = "#F18204",
"wk40" = "#B56203"),
soy = c("wk0" = "#FE71B1", "wk10" = "#FE318E", "wk20" = "#F4016E",
"wk30" = "#B70153", "wk40" = "#7A0137")
)
create_crop_line_plot <- function(data, crop_name, color_palette) {
# Filter data for specific crop
crop_data <- data %>%
group_by(timepoint) %>%
filter(crop == crop_name)
# Ensure absorbance is numeric
crop_data$normalized_absorbance <- as.numeric(crop_data$normalized_absorbance)
# Create the plot
# Scale the data to match ridge plot proportions
max_abs <- max(crop_data$normalized_absorbance, na.rm = TRUE)
line_plot <- crop_data %>%
ggplot(mapping = aes(x = wavenumber,
y = normalized_absorbance,
group = filename,
color = timepoint)) +
xlim(3200, 1200) +
geom_line(alpha = 0.8) +
theme_classic() +
scale_color_manual(values = color_palette) +
ggtitle(paste(str_to_title(crop_name), "root suberin")) +
labs(x = "Wavenumber (cm⁻¹)",
y = "normalized absorbance") +
# Suberin peak annotations
# Aliphatic suberin region (3000-2800 cm-1)
annotate("rect", xmin = 2800, xmax = 3000, ymin = -Inf, ymax = Inf,
alpha = 0.1, fill = "#5076ab") +
# Specific aliphatic peaks
geom_vline(xintercept = 2976, linetype = "longdash", linewidth = 0.3, color = "#3D6AA9") +
geom_vline(xintercept = 2959, linetype = "longdash", linewidth = 0.3, color = "#3D6AA9") +
geom_vline(xintercept = 2930, linetype = "longdash", linewidth = 0.3, color = "#3D6AA9") +
geom_vline(xintercept = 2856, linetype = "longdash", linewidth = 0.3, color = "#3D6AA9") +
# Suberin esters (1737-1740 cm-1)
annotate("rect", xmin = 1733, xmax = 1741, ymin = -Inf, ymax = Inf,
alpha = 0.15, fill = "#4d653a") +
geom_vline(xintercept = 1737, linetype = "dashed", linewidth = 0.3, color = "#4d653a") +
# Aromatic suberin (1506 cm-1)
geom_vline(xintercept = 1506, linetype = "dotted", linewidth = 0.3, color = "#5f3e38") +
# Add aliphatic annotation
annotate("text", x = 2600, y = max(crop_data$normalized_absorbance)-0.05, label = "aliphatic suberin",
size = 3, hjust = 0.5) +
# Add esters annotation
annotate("text", x = 1910, y = max(crop_data$normalized_absorbance), label = "suberin esters",
size = 3, hjust = 0.5) +
# Add aromatic annotation
annotate("text", x = 1300, y = max(crop_data$normalized_absorbance)+0.02, label = "aromatic suberin",
size = 3, hjust = 0.5)
return(list(plot = line_plot, data = crop_data))
}
# Generate plots and tables for each crop
crops <- c("noPlant", "wheat", "rice", "soy")
for (crop in crops) {
# Create plot and get data
result <- create_crop_line_plot(dpt_data_roots, crop, crop_timepoint_colors[[crop]])
# Print the plot
print(result$plot)
cat("<br>") # Add space between plot and table
# Create and print ID combinations table
id_table <- create_crop_id_table(result$data, crop)
cat(paste0("<h4>", str_to_title(crop), " pot root samples</h4>"))
cat(knitr::kable(id_table, format = "html", escape = FALSE, table.attr = 'style="border-collapse: separate; border-spacing: 35px 0;"'))
cat("<br><br>") # Add HTML line breaks for spacing
}| wk30 | wk40 |
|---|---|
| 052 | 052 |
| wk0 | wk10 | wk40 |
|---|---|---|
| 001, 032 | 029, 053 | 029, 053, 109 |
| wk0 | wk10 | wk40 |
|---|---|---|
| 103, 107 | 055, 087, 119 | 055, 087, 119 |
| wk0 | wk10 | wk40 |
|---|---|---|
| 050, 062 | 012, 080, 094 | 030, 080, 094 |
## Summary of FTIR Analysis:
## Number of samples: 32
## Crops: wheat, soy, rice, noPlant
## Timepoints: wk0, wk10, wk40, wk30
## [1] "Crop Statistics:"
## # A tibble: 4 × 7
## crop Mean_simple_plant_prop SE_simple_plant_prop Mean_complex_plant_prop
## <fct> <dbl> <dbl> <dbl>
## 1 noPlant 0.397 0.0354 0.205
## 2 wheat 0.430 0.0103 0.192
## 3 rice 0.451 0.00973 0.182
## 4 soy 0.368 0.0292 0.199
## # ℹ 3 more variables: SE_complex_plant_prop <dbl>, Mean_microbial_prop <dbl>,
## # SE_microbial_prop <dbl>
## [1] "\nTimepoint Statistics:"
## # A tibble: 4 × 7
## timepoint Mean_simple_plant_prop SE_simple_plant_prop Mean_complex_plant_prop
## <chr> <dbl> <dbl> <dbl>
## 1 wk0 0.436 0.0204 0.181
## 2 wk10 0.398 0.0313 0.190
## 3 wk30 0.362 NA 0.216
## 4 wk40 0.412 0.00620 0.203
## # ℹ 3 more variables: SE_complex_plant_prop <dbl>, Mean_microbial_prop <dbl>,
## # SE_microbial_prop <dbl>